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CONTRACTOR REPORT 


COMPUTATION OF TURBULENT FLOWS USING AN EXTENDED 
k-e TURBULENCE CLOSURE MODEL 


INTRODUCTION 


It is well known that, in the framework of isotropic or eddy viscosity turbulence 
closure modeling, good predictions of the mean and turbulent flow fields rely on 
reasonable descriptions of the turbulent length scale and velocity scale inside the flow 
field. The most widely used isotropic two-equation turbulence model is the standard 
k-e model of Reference 1. This model, with standard model constants, has been 
applied over a wide range of turbulent flow problems [2,3,4]. In general, the stand- 
ard k-e model gives qualitative predictions for most turbulent flows, while the well- 
known plane jet and round jet anomaly is one of the typical problems that the standard 
k-e model had failed to give consistent predictions as was described by Pope [5] and 
Hanjalie and Launder [6], Although the spreading rate of plane jet was predicted 
correctly by the standard k-e model, the spreading rate of round jet was over pre- 
dicted by more than 30 percent. Also, for some bench-mark elliptic flow problems, 
such as a turbulent flow over a backward- facing step [4,7] and a turbulent confined 
swirling flow problem [8], etc., the standard k-e model gives highly diffusive results. 
The reattachment length of the recirculation region of the backward- facing step flow, 
for instance, was measured to be around 7.2 step heights. However, the consensus 
prediction of the reattachment length using the standard k- e model was reported to 
be around 5.2 step heights [7]. 

The inconsistency of the standard k-e model is very often attributed to the 
dissipation rate equation which is highly empirical in nature. Improvement of the 
model performance is usually achieved by modifying the dissipation rate equation. In 
References 5 and 6, extra irrotational strains are added to the production term of the 
dissipation rate equation. With these modifications, some improvements in the generality 
of the k-e model have been demonstrated [9,10]. For internal recirculating flow 
problems, many implementations for the k-e model have been proposed by employing 
the so-called Richardson number correction for the e-equation [11,12], Although some 
improved predictions for certain types of elliptic flow problems have been reported 
using the Richardson number corrections, they were found to be problem dependent 
in nature [11,12]. 

In the present paper, a similar but more general approach is taken by adding 
a second time scale of the production range of turbulence kinetic energy spectrum to 
the dissipation rate equation. This extra time scale enables the energy transfer 
mechanism of the turbulence model to respond to the mean strain more effectively. 

The present modification results in one extra term along with one extra modeling con- 
stant added to the standard k-e model. The new set of model constants are tuned 
based on the experimental data of homogeneous turbulence decay with or without mean 
strain imposed [13,14,15], a near wall boundary conditions analysis and numerical 
optimization. 

Numerical examples included in the present paper involve fully developed tur- 
bulent channel and pipe flows [16,17], turbulent free-shear flows [18,19,20,21], flat 
plate turbulent boundary layer [22], turbulent flow over a backward- facing step [4,7] 



and a confined turbulent swirling flow [8], This wide variety of example problems is 
included to assess the performance and generality of the present extended k- e tur- 
bulence model. For numerical computations, a second-order accurate finite difference 
boundary layer code was used for boundary layer type flow problems and a nearly 
second-order accurate finite difference elliptic flow solver was used for the internal 
turbulent flow problems in the present study. The boundary layer code employs t 
stream- wise marching technique using a second-order backward differencing scheme. 
In the transverse direction, a second-order upwind differencing scheme is used for 
the convection terms and a central differencing scheme is used for the diffusion 
terms. In the elliptic flow solver, the diffusion terms of the governing equations are 
discretized by using a central differencing scheme and the convection terms are dis- 
cretized by using a second-order upwind differencing scheme [23,24]. Favorable 
results of many other example problems using the present k- e turbulence model and a 
rigorous and accurate finite element boundary layer code can be found in Reference 
25. Results of Reference 25 reflect and confirm the consistency of the present tur- 
bulence model using both finite difference and finite element methods. 


GOVERNING EQUATIONS 


For steady incompressible turbulent flows, the flow field can be characterized 
by conservation laws. These are the continuity equation and the Navier- Stokes (or 
momentum) equations. With the commonly used Reynolds decomposition and averaging 
procedure, the system of governing equations for turbulent flow field can be written 
in the following form [ 26] . 


(pu) x + (pv) y = 0 


( 1 ) 


PU u x tpTUy- (u e V x ■ (ti e u y ) y = _p x + ( Vx u x + ^eVx 


PU V + pv V^ ~ (y V ) - (u V ) = -p + ( u ) U + (u ) V 

p x p Vt ^e x'x y y *y M e'x y e y y 


(2a) 

(2b) 


y e = » + p t 


where u and v are velocity vectors in x- and y-coordinate directions, respectively, 
y and p represent the fluid viscosity and density, respectively, and p denotes the 
static pressure. The above governing equations can be simplified for boundary layer 
type flows [27]. In the above equations, the Boussinesq assumption has been used 
to relate the Reynolds stresses to the mean strains through the turbulent viscosity, 
y^. The turbulent eddy viscosity is then related to the turbulence length scale, l^, 

and the turbulence velocity scale, u^., which can be expressed in terms of other tur- 
bulence quantities, namely the turbulent kinetic energy, k, and its dissipation rate, 
e. That is. 


n t = p l t u t 



k 3/2 /e 


u. 


= k 1/2 


(2c) 


2 



where C is usually taken as a constant, = 0.09. Although previous investiga- 
tions [28,29,30] have suggested that should be a function of flow field quantities, 

modification of C will not be considered here. Effects of C function on the per- 

y y 

formance of turbulence model can be found in References 25 and 30. The turbulence 
quantities, k and e, are provided by employing an extended k-e turbulence closure 
model which is described in the following section. 


EXTENDED k-e CLOSURE 


For the closure of the governing equations, the widely used two-equation k-e 
turbulence modeling approach is employed in the present study. The transport equa- 
tion of the turbulent kinetic energy can be written as [ 1] : 

pu k x + pv k y - [(y t /a k ) k x ] x - [(y t /c k > k y ] y = p(P r - e) (3) 

where and e represent the production rate and the dissipation rate of the turbu- 
lent kinetic energy, k, respectively, and o k is a modeling constant. The production 

rate is related to the mean strain of the velocity field through the Boussinesq assump- 
tion . That is , 

P r = (li t /p)[(u y + v x ) 2 + 2(u x 2 + v y 2 )] . 


For boundary layer type flows, the production rate can be simplified to be: 
P r = (Pj/P) U y 2 • 


However, for rapidly evolving flow fields (e.g., plane jet or round jet exhausting 
into a still air) , it is more appropriate to retain the full production term to have 
better representation of the energy generation rate. 

For the dissipation rate equation, two time scales are included to allow the dis- 
sipation rate to respond to the mean strain more effectively than that of the standard 
k-e model. This is the major improvement of the present k-e model for complex tur- 
bulent flow problems. The time scales included in the present model are: the pro- 

duction range time scale, k/P^, and the dissipation rate time scale, k/e. The final 

expression of the dissipation rate transport equation is given as: 


PU £ x + pv £y - [(U t /O e ) C x ] x - t <U t /O e ) Ey ]y = 

p(C l P r e/k ■ c 2 e2/k + C 3P r 2/k) • 


(4) 
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The last term of equation (4) represents the energy transfer rate from large scale 
turbulence to small scale turbulence controlled by the production range time scale and 
the dissipation rate time scale. The net effect of the present energy transfer function 
enhances the development of e when the mean strain is strong, or large production 
rate, and the generation of e is suppressed when the mean strain is weak, or small 
production rate. Consequently, as the model constants are carefully tuned, the * 
present formulation enables the dissipation rate to respond to the mean flow field more 
rapidly so as to control the development of the turbulent kinetic energy more effec- 
tively . 


To determine the model constants, C^, C 2 , and C^, experimental data of the 

decay of homogeneous turbulent flows, with or without mean strains [13,14,15], and 
a simplified near wall analysis [25] are used to decide the feasible range of these 
constants. Numerical optimization is then employed to anchor the final model constants 
by matching the predictions to the measured data of several simple turbulent flows, 
i.e., fully developed turbulent channel and pipe flows and turbulent plane jet. For 
the rest of the two model constants, a^. and a £ , two simple physical arguments are 

imposed to set up two criteria for determining these two constants. First, any tur- 
bulence model must satisfy the realizability [31], especially near the edge of boundary 
layers, such that the dissipation rate would be vanishing, in the direction away from 
the boundary layer, at a rate faster than that of the turbulent kinetic energy. This 
requires that o k is smaller than a . Secondly, by observing many data of boundary 

type turbulent flows [18,19,20,21,22], it is obvious that the boundary edges of the 
turbulent kinetic energy extend much wider (about 10 percent wider) than that of the 
mean velocity profiles. This implies that it is proper to have a, less than unity. 

The final numbers of these model constants are: K 


a k = 0.75 


a £ = 1.15 ; C 1 = 1.15 ; C 2 =1.9 ; C 3 = 0.25 . 


Detailed descriptions about the determination of these model constants can be found 
in Reference 25. 


NUMERICAL METHOD AND BOUNDARY CONDITIONS 


A finite difference marching procedure is employed in the present study to 
carry out the computations of boundary layer type turbulent flow problems. The 
present finite difference method includes a second order central difference approxima- 
tion for the cross-stream diffusion terms. For the cross-stream convection terms, a 
second order upwind differencing scheme is employed. In the stream- wise (or march- 
ing) direction, a second order backward (or upwind) differencing scheme is used to 
provide a fully implicit boundary layer marching procedure. The whole numerical 
scheme has the accuracy of second order. Uniform grid systems in the longitudinal 
and transverse directions are employed to cover the entire physical domain. This 
requires that the initial computational domain in the transverse direction extend 
outward several boundary layer thicknesses. 

At the initial data station, experimental measurements of the mean velocity pro- 
file and the turbulent kinetic energy distribution are employed. Measurements of 
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Reynolds stress are used to generate the initial dissipation rate distribution through 
the Boussinesq assumption and equation (2c). Along the free stream boundary and 
the center line of jet flows, zero normal gradients (or vanishing flux) boundary con- 
ditions are imposed. Near the solid wall boundary, the commonly used wall function 
approach [2] is employed to provide boundary conditions for the velocity and tur- 
bulence quantities. This wall function approach is based on the assumptions of a* 
logarithmic wall law velocity profile and a near wall turbulent kinetic energy equilib- 
rium condition. Although these assumptions are not valid for boundary layers subject 
to stream- wise curvature effects, strong adverse pressure gradient or flow separation, 
the choice of the approach for solid wall boundary conditions is, to the best knowledge 
of the authors, due to the lack of a better and more general alternative at the pre- 
sent time. 

For elliptic flow problems, a 2- dimensional finite difference elliptic flow solver 
is employed [23,24]. In this elliptic flow solver, central differencing is used to 
approximate the diffusion terms. For the convection terms, a second order upwind 
differencing scheme is implemented in the basic solver of References 23 and 24 to 
provide numerical stability and retain numerical accuracy close to second order. The 
pressure field is obtained by employing the SIMPLE-C algorithm of Reference 32, 
which corrects the pressure and velocity fields successively until the conservation of 
mass is satisfied throughout the entire computational domain. The convergence cri- 
terion of the present elliptic flow solver requires the sum of the maximum correction 
of the velocity and pressure between two successive iterations to be less than 10" 4. 

At the entrance of the computational domain, inlet boundary conditions are 
specified. At the exit, zero gradient boundary conditions are imposed. Near the 
solid wall, the conventional wall function approach is used [2]. For axisymmetric 
flow problems, symmetric boundary conditions are imposed along the axis. 


RESULTS AND DISCUSSION 


In order to assess the performance of the present turbulence model, several 
boundary layer and elliptic turbulent flow problems were tested in the present study . 
These flow problems include: fully developed turbulent channel and pipe flows [16,17]; 

turbulent plane jet and round jet [18,19,20,21]; a flat plate turbulent boundary layer 
flow [22]; an internal turbulent recirculating flow over a backward facing step [4,7]; 
and a confined turbulent swirling flow [8]. Results of the computation of the above 
flow problems are discussed in detail below. 

1. Fully Developed Channel and Pipe Flows 

For these test cases, experimental data given by Laufer [16,17] were used for 
data comparisons. The Reynolds numbers of the selected cases of the fully developed 
channel and pipe flows were 61,600 and 200,000, respectively. The center line velocity 
of the channel and pipe flows were 7.05 m per sec and 30.05 m per sec, respectively. 
For fully developed flow problems, the boundary layer equations can be further 
simplified to be a system of ordinary differential equations. Symmetric boundary con- 
ditions were specified at the channel or pipe center line. Near the wall boundary, 
experimental data of the mean velocity and the turbulent kinetic energy were speci- 
fied as fixed boundary conditions which were intended to avoid errors arising from 
the use of wall function approximations. The near wall dissipation rate was estimated 
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by using the mixing length approximation. The longitudinal pressure gradient was 
given by the experimental data. Fifty grids were used for the computations of these 
cases . 


Figures 1(a) and 1(b) show comparisons of the predicted and the measured 
mean velocity and turbulence quantities for the channel and pipe flows, respectively. 

It is clear in Figure 1 that the present turbulence model gives satisfactory predictions 
for these test cases. The standard k-e model gives almost the same predictions. 

2. Submerged Jets 

Two submerged- jet , i.e., jet exhausting into a still air environment, problems 
were studied which included a submerged plane jet [18] and a submerged round jet 
[19] . Spreading rate of the jet half width, dy^^/dx, is of major concern here. Experi- 
mentally, the submerged jet can reach a constant spreading rate within 30d downstream 
of the jet exit, where d stands for the jet nozzle width or diameter. The experimental 
spreading rates of the submerged plane jet and round jet were found to be 0.11 and 
0.098, respectively. 

In the present study, computations of the submerged jets started from lOd 
downstream of the jet exit. The initial velocity and turbulent kinetic energy profiles 
were generated using the similarity profiles. The initial data for the dissipation rate 
were then estimated from the Reynolds stress distribution through the Boussinesq 
assumption and equation (2c). Computational domain in the transverse direction 
extended outward 10 times the initial jet width to allow marching to 40d downstream. 
Due to the rapidly evolving velocity field, the full production expression was used in 
the present turbulence model to give a better representation of the non- boundary- 
layer character of the flow field. This kind of implementation was found to be only 
necessary for the submerged jet problems since velocity gradients in the stream- wise 
direction are not negligible. 

The spreading rates of the submerged jets were then calculated from the solu- 
tions of the marching procedure. The results show that the computed spreading 
rates of the submerged plane and round jets are 0.111 and 0.108, respectively. 

Using the same form of production term, the standard k-e turbulence model produced 
spreading rates of 0.117 and 0.119 for the plane and round jets, respectively. With 
the simplified production term (as commonly used), the standard k-e model gave 
spreading rates of 0.11 and 0.125 for the two jets, respectively. From these results, 
it is apparent that the effect of the full production term treatment is significant for 
rapidly developing flow fields. 

3. Jets Exhausting into Moving Stream 

Two co-flowing jet problems were included in the present study. These two 
cases correspond to the experimental investigations of Bradbury [20], for a plane jet 
exhausting into a moving stream of speed 7.344 m/s, and of Antonia and Bilger [21], 
for an axisymmetric jet exhausting into a moving stream of speed 30.5 m/s. For the 
plane jet problem, the jet velocity at the exit of a jet slot with width, d = 0.009525 m, 
was 45.9 m/s. Numerical computation was started from lOd downstream of the jet 
exit, where measured data of mean velocity and turbulence intensities distributions 
were given. Solutions were obtained by using the present marching procedure from 
x = lOd to x = 70 d. Ninety-one grids in the transverse direction were used in the 
marching procedure. Figure 2 illustrates comparisons of the measured and the 
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predicted decay of the center line velocity, u c , and growth of the jet half width, 

y 1 ^ 2 * Detailed comparisons of the mean velocity profiles and the distributions of 

turbulence quantities are shown in Figure 3 for x = 40d and x = 68d. It is clear from 
these results that the present turbulence model reproduces the experimental data very 
well except slight discrepancies in the turbulent kinetic energy profile at x = 68d^ 

For the case of round jet, the jet exit velocity out of a nozzle with diameter, 
d = 0.00528 m, was 137 m/s. The same marching procedure was used for this problem. 
Since only self-similar profile data were given in Reference 21, the initial data were 
generated from the similarity profiles which do not resemble the real test conditions. 
Therefore, only the comparisons of the self-similar profiles are presented in Figure 4. 
Numerically, self-similar profiles were produced using results after 50d downstream of 
the initial data plane. Figure 4 shows that the present turbulence model also gives 
good predictions for this round jet problem. 

4. Flat Plate Turbulent Boundary Layer 


In this example, a turbulent boundary layer developing along a flat plate with 
zero pressure gradient was examined. The experimental case of Wieghardt [22] was 
simulated. The free stream velocity was 33 m/s. The initial data plane was located 
at x = 0.94 m and the marching procedure was performed up to x = 5.0 m. Ninety- 
one grids in the transverse direction were used in the marching procedure. Figure 5 
shows the comparisons of wall shearing stress distributions using the standard k- e 
model and the present extended k-e model. It can be seen from Figure 5 that the 
standard k- e model over-predicts the wall shearing stress by about 5 percent while 
the present turbulence model under-predicts the wall shearing stress by about the 
same percentage. Figure 6 illustrates the comparisons of mean velocity profiles, 
turbulent kinetic energy profiles and Reynolds stress distributions at two stations, 
i.e., x = 2.887 m and x = 4.987 m. It is clear in Figure 6 that the mean velocity 
profiles are well predicted by the present turbulence model while discrepancies in 
turbulence quantities are apparent. This is mainly attributed to the inadequacy 
of the wall boundary conditions for k and c using the conventional wall function 
approach [2] which gives low turbulent kinetic energy at the wall function point. 
Similar results were obtained using the standard k-e turbulence model. Better wall 
boundary conditions may improve the predictions of the turbulence quantities in both 
models. 

5. Backward- Facing Step Flow 

One of the standard test cases of complex elliptic turbulent flow presented in 
the Stanford Conference [4] was a confined turbulent recirculating flow over a 
backward- facing step [7]. The experimental configuration contained a straight 
channel followed by a sudden expansion with 2:3 expansion ratio. The inlet flow 
velocity was about 18 m/s. The measured size of the separation region (reattachment 
length) caused by the sudden expansion was around 7.2 step heights. The computa- 
tion of the reattachment length of this problem has been used extensively for assess- 
ing the performance of turbulence models [4], 

In the present study, the computational domain extended 6 step heights upstream 
of the expansion plane and extended 30 step heights downstream of the expansion 
plane to enable the application of zero gradient exit boundary conditions. At the 
flow inlet, uniform velocity and turbulent kinetic energy profiles were specified. 

Inlet boundary condition of the dissipation rate was then estimated by assuming a 
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constant mixing length equivalent to 0.03 times the inlet channel width. Near solid 
wall boundaries, the conventional wall function approach was employed to provide 
boundary conditions for \the momentum equations and the turbulence model. A 51 by 
41 grid system was used for this case. This grid size yields grid independent solu- 
tions in which the difference in the reattachment length using grid sizes of 41 by 35 
and 51 by 41 was less than 1 percent. The standard k-e model and the present k- e 
model were used for the current test case. 

Figure 7 shows the comparisons of the locus of flow reversal (where the longi- 
tudinal velocity component, u, changes sign). It is clear that the present turbulence 
model gives much better predictions than the standard k-e model. The present tur- 
bulence model predicts a reattachment length of 7.0 step heights while the standard 
k-e model gives only 5.2 step heights. Figure 8 illustrates the comparisons of wall 
static pressure distributions along the step side wall. Good agreement between the 
present model predictions and the measured data is shown clearly in Figure 8 while 
apparent discrepancies are revealed for the standard k-e model predictions. Figures 
9(a) and 9 (b) show detailed comparisons of the mean velocity profiles and the dis- 
tributions of the turbulent kinetic energy downstream of the step. Again, improve- 
ments of the present k-e model over the standard one are illustrated in the region 
near reattachment. It can be seen from Figure 9(b) that the turbulent kinetic energy 
overshoot phenomenon downstream of the expansion plane of the standard k-e model 
is suppressed effectively by the present turbulence model which contains a more 
effective energy transfer function in the e-equation, equation (4), see x/h = 2.33. 

6. Confined Swirling Flow 

In the present study, a confined turbulent swirling flow problem investigated 
experimentally by Roback and Johnson [8] was solved numeric ally. Several numerical 
studies of this flow using the standard k-e turbulence model have been reported 
[32,33, 34] . Experimentally , the swirling flow channel consisted of two coaxial inlet 
pipes with the swirling guide vanes installed between the inner and outer pipes. The 
inner and outer pipes had radii of 12.5 mm and 29.5 mm, respectively. The inlet 
channel was followed by a sudden pipe expansion with expansion ratio around 1:2. 

The radius of the downstream pipe, R q , was 61 mm. The inlet swirling velocity 

generated by the swirling guide vane created a central recirculation zone along the 
pipe center line downstream of the expansion plane. This central recirculation zone 
was accompanied by a corner recirculation region downstream of the step. 

In order to avoid geometrical and flow complexities upstream of the expansion 
plane, where a three-dimensional flow field was expected and detailed experimental 
data were not available, the present numerical computation employed the first data 
plane, which is 5 mm downstream of the expansion plane, as the inlet boundary where 
detailed experimental data were provided [8]. The exit boundary was located 1800 mm 
downstream of the expansion plane. A grid size of 51 by 41 was used for this test 
case. More grids were clustered near the inlet where rapid flow development is 
expected. The conventional wall function approach was again used for near wall 
boundary conditions. A converged solution was obtained in 1500 iterations due to the 
use of a secohd order upwind differencing scheme for the convection terms. 

Computed results of the sizes of the central recirculation zones from the use of 
the present k-e model and the standard k-e model are compared with the measured 
data illustrated in Figure 10. Comparisons of the development of the center line 
velocity are presented in Figure 11. It is clear in Figures 10 and 11 that the present 
turbulence model predicts much better results than that of the standard k-e model. 
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Figure 12 gives a detailed comparison of the radial distribution of the axial velocity 
at x = 25 mm. Again, the present turbulence model represents the measured data 
much better. Also, results predicted by the standard k-e model are similar to those 
reported previously in References 32, 33, and 34. 


CONCLUSIONS 


In the present paper, an extended k-e turbulence model has been proposed and 
tested for a wide variety of turbulent flow problems. An extra time scale has been 
introduced in the construction of the transport equation of the rate of dissipation of 
the turbulent kinetic energy. This results in an energy transfer function which has 
been found to be more effective than that of the standard k-e turbulence model. 
Numerical tests of example problems including simple boundary type turbulent flows 
and complex internal turbulent recirculating flows have demonstrated that the present 
model is more general than the standard k-e model. For simple boundary type 
turbulent flows, the present model gives similar results predicted by the standard k-e 
model. However, for complex elliptic turbulent flow problems, which involve rapid 
changes of turbulent kinetic energy production and dissipation rates, the present 
extended k- e has been shown to give much better results than the standard k- e 
model. This is mainly attributed to the energy transfer function introduced in the 
dissipation rate equation which enables the development of the field of the dissipation 
rate to suppress the overshoot phenomenon of the turbulent kinetic energy , as 
revealed by the standard k-e model, when the mean shear is strong. 

Applications of the present model to other complex turbulent flow problems, 
e.g., turbulent flows inside strongly curved ducts, turbulent flows inside turbo- 
machinery passages, time dependent turbulent flow problems, etc., are planned for 
future study. This will give a better assessment of the generality of present tur- 
bulence model. 
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Figure 1. Fully developed channel and pipe turbulent flows. 
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Figure 5. Wall shearing stress along a flat plate 
with zero pressure gradient. 
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Figure 6. Flat plate turbulent boundary layer. (See Fig. 1 for legend.) 
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Figure 7. Locus of flow reversal downstream of a backward- facing step. 
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Figure 8. Static pressure coefficient distribution along the step side wall 
of a backward -facing step turbulent flow. 
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Figure 9. Mean velocity and turbulent kinetic energy profiles 
downstream of a backward-facing step. 
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Figure 10. Envelope of the central recirculation zone 
of a confined swirling turbulent flow. 



Figure 11. Center line axial velocity distribution of a 
confined swirling turbulent flow. 
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Figure 12. Radial distribution of the axial velocity at x = 25 mm 
downstream of the coaxial swirling jets. 
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